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The Born effective charge, Z* , that describes the polarization created by collective atomic dis- 
placements, can be computed from first-principles following different techniques. We establish the 
connections existing between these different formulations, and analyze the related band-by-band de- 
compositions. We show that unlike for the full Z* , the different band-by-band values are not equal, 
and emphasize that one of them has a natural physical meaning in terms of Wannier functions. 



I. INTRODUCTION 



The Born effective charge (Z*) is the equivalent, for crystalline insulating solid, of the atomic polar tensor 
introduced for isolated molecules . It is a dynamical charge related to the macroscopic polarization induced by the 
collective displacements of nuclei belonging to a given sublattice. In the study of the lattice dynamics of insulating 
crystals, it is considered as a fundamental quantity, because it governs the amplitude of the long-range Coulomb 
interaction between nuclei, and the splitting between longitudinal (LO) and transverse (TO) optic phonon modes. 

In simple materials, like B^~^ binary crystals the phonon eigenvectors are imposed by symmetry. Infra- 
red measurements of the splitting between LO and TO modes allows an accurate estimation of \Z*\'^ /too and offers 
therefore an unambiguous way to extract the amplitude of Z* from the experiment. However, in more complex 
materials like ABO3 compounds, LO and TO mode eigenvectors are not necessarily equivalent. The determination of 
Z* from the experimental data is consequently not straightforward and requires the use of some approximations [Q. 
For such compounds, the development of theoretical methods giving direct access to Z* acquires therefore a specific 
interest. 

Conventionally, the Born effective charge tensor Z* of nuclei belonging to the sublattice k is defined as the 
coefficient of proportionality relating, under the condition of zero macroscopic electric field, the change in macroscopic 
polarization Vfi along the direction [3 and the collective nuclear displacements of atoms k along direction a, times the 
unit cell volume f2o : 
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However, a thermodynamical equality relates the macroscopic polarization to a derivative of the electric enthalpy E 
and another relationship connects the forces Fi^ on the nuclei k to a derivative of the electric enthalpy so that Z* can 
be alternatively defined as follows: 
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From these relationships, Z* can therefore be thought either as (i) the change of polarization induced by the collective 
displacements of atoms k, under zero field (ii) a mixed second derivative of the electric enthalpy or (iii) the derivative 
of the force induced on a nucleus k, by an homogeneous effective electric field Efj, at zero atomic displacements. 

The three definitions are formally equivalent. However, each of them can lead to different algorithms for the 
computation of Z* from first-principles. Among the approaches which are the most widely used, a first powerful and 
systematic procedure was introduced by Baroni, Giannozzi and Testa |^, who suggested to determine Z* from a linear 
response formalism grounded on a Sternheimer equation. A different algorithm, based on a variational principle, was 
later reported by Gonze, Allan and Teter Q, yielding a new alternative expression for Z* . Thanks to progress in the 
theory of the macroscopic polarization, Z* is also directly accessible from finite difference of polarization The first 
two algorithms were exclusively implemented within the density functional formalism (DFT) while the last one also 
allowed calculations of changes in polarization within different other one-electron schemes (Hartrec-Fock method [|j , 
model GW approximations to many-body theory [oljiol], Harrison tight- binding model jllj) and the Hubbard tight- 
binding model 
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Accurate predictions of the Born effective charges have been reported for a large variety of materials. In some of 
these studies |9|,p"3|-pT[ , the decomposition of Z* in individual contributions from separate groups of occupied bands 
appeared as a powerful tool to identify the microscopic mechanisms monitoring its amplitude. However, the physical 
interpretation of these decompositions was never explicitely discussed. Moreover, unlike the total Z*, contributions 
from individual groups of bands are not uniquely defined. 

In the present paper, we aim at presenting the links between the theoretical frameworks used nowadays to compute 
Z*, and at deducing from this comparison the correct way to develop a band- by-band analysis. We demonstrate 
that the natural decompositions arising from Eq. (|l|) and Eq. dijfers, although independent approches (linear 
response. Berry phase approach, electronic Wannier functions) to Eq. (|^) are strictly equivalent provided the phase of 
the wavefunctions are correctly chosen. 

The paper is organized as follows. In Section II, we describe the mathematical links existing between the different 
expressions that can be used to determine the global value of Zi^^ap within the density functional formalism without 
yet referring to band- by-band decompositions. We adopt the notations of Ref. In Section III, we discuss how 

contributions from isolated sets of bands can be separated from each others. We identify different expressions and 
discuss their meaning in terms of electronic Wannier functions. In Section IV, we illustrate our results on a numerical 
example, emphasizing that independent decompositions yield in practice radically different values. Finally, in Section 
V, we conclude on the physically correct way of performing band- by-band decomposition of Z*. 



II. DIFFERENT FORMULATIONS OF Z* 



For practical purposes, in what follows, the Born effective charge will be decomposed into two contributions : 

^K.aP — ZkSoP + Z^^^p. (3) 

The first term, Z^, is the charge of the nuclei (or pseudo-ion, in case of pseudopotential calculations), and can be 
trivially assigned. The second, Z'^ is the contribution due to the electrons. 



A. First derivative of the polarization 

A first straightforward approach for the determination of Z* consists in computing the difference of macroscopic 
polarization between a reference state, and a state where the atoms belonging to the sublattice n have been displaced 
by a small but finite distance Atk,q. The electronic contribution to Z* can be obtained as : 

In periodic systems, the change in electronic polarization in zero field can be computed from the King-Smith and 
Vanderbilt formula : 

-j^ OCC n Q 

where s is the occupation number of states in the valence bands (s = 2 in spin-degenerate system) and u„k is the 
periodic part of the Bloch functions. Taken independently, the matrix elements of the previous equation are ill-defined 
because the phase of the wavefunctions at a given wavevector of the Brillouin zone is arbitrary, and thus unrelated 
with the phases at neighbouring k points. However, the integral of the right-hand side is a well-defined quantity, 
which takes the form of a Berry phase of band n, as discussed by Zak | pO[ |. 

The King-Smith and Vanderbilt definition is valid only under the constraint that the wavefunctions fulfill the 
periodic gauge condition. This means that the periodic part of Bloch functions must satisfy 

Unk(r) = e*'^-'' w„k+G(r), (6) 

This condition does not fix unambiguously the phase of the wavefunctions at a given k-point (even not at neighbouring 
k-points) but it imposes a constraint for wavefunctions at distant wavevectors. It defines a topology in k-space, within 
which the polarization takes the convenient form of a Berry phase. 
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When working within one-electron schemes (DFT, Hartree-Fock, ...), a second choice of phase is present at another 
leveL For the ground-state, the Lagrange multipUer method appUed to the minimization of the Hohenberg and Kohn 
fonctional under orthonormahzation conditions on the wavefunctions , gives the following equations : 

occ 

-ffklWmk) = ^ A„„,k \Unk) (7) 

n 

This condition, associated with the minimisation of the Hohenberg and Kohn energy functional, means that the 
application of the Hamiltonian to a given wavefunction generates a vector which must stay within the Hilbert space 
defined by the set of Unk wavefunctions. We observe that a unitary transformation between the wavefunctions will 
leave that Hilbert space invariant, and Eq. ^ will remain satisfied provided the matrix of Lagrange multiplier Amn.k 
is transformed accordingly. In order to build Kohn-Sham band structures, the unitary transform is implicitely chosen 
such as to guarantee 

Amn,k — ^mn ^m,k (8) 

in which case 6j7i.k correspond to the eigenvalues of the Kohn-Sham Hamiltonian and the associated functions iid.miL 
are the Kohn-Sham orbitals. This choice is called the diagonal gauge condition. All along this work, it will be 
emphasized by a "d" subscript. 

We note that the periodic gauge condition connect wavefunctions at different k-points, while the diagonal gauge 
condition fixes wavefunctions at a given k-point. The choice defined by Eq. is not mandatory, and the computation 
of the total energy, the density, or the Berry phase (Eq. (||)) will give the same value independently of the fulfillment 
of Eq. (H). The diagonal gauge is the natural choice for the ground-state wavefunctions while, as it will be discussed 
later, another choice is usually preferred for the change in wavefunctions in linear-response calculations. 

Instead of approximating Eq. (|l|) from finite differences, it can be computed directly. The combination of Eqs. (|^), 
(I) and §) gives : 

el _ [ r,9u„k|9u„k, , , I 9 ,(9u„k, 



where the second expectation value can be worked out : 

In the last expression, the first term of the right-hand side is the gradient of a periodic quantity integrated over the 
Brillouin zone. Within any periodic gauge, its contribution will be zero. Using the time-reversal symmetry, we arrive 
therefore at the final expression : 

el _ r, . f , <9u„k I 9u„k 
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The first-derivatives of the wave functions, 9unk/9T„_ct and dunk/dkp, appearing in this expression, can be computed 
by linear-response techniques either by solving a first-order Sternheimer equation P,|2l| or by the direct minimization 
of a variational expression as described in Ref . [^,0 . 

We note that the choice of gauge will influence the value of the first-derivative of it„k, although the integrated 
quantity Z'^^p must remain independent of this choice (in any periodic gauge). Usually, the following choice is 
preferred in linear-response calculations : 



|Umk) = (12) 



P 



for m and n labelling occupied states, and A representing either the derivative with respect to the wavevector or to 
atomic displacements. As emphasized by the "p" subscript, this condition defines what is called the parallel gauge 
and insures that the changes in the occupied wavefunctions are orthogonal to the space of the ground-state occupied 
wavefunctions. This projection on the conduction bands is not reproduced within the diagonal gauge defined by the 
generalization of Eq. (0]|) at the first order of perturbation, as elaborated in Ref. [p2| . 
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B. Mixed second derivatrive of the electric enthalpy 



The Born effective charge also appears as a mixed second derivative of the electric enthalpy. Therefore, as reported 
in Eq. (41) of Ref. ||l9|, Z^;} can be alternatively formulated in terms of a stationary expression, involving the 
first-order derivative of the wavefunctions with respect to a collective displacement of atoms of the sublattice k, and 
the first-order derivatives of the wavefunctions with respect to an electric field and to their wavevector [p3| . 

The mathematical equivalence between Eq. (^ij) and Eq. (41) of Ref. |l^ is a consequence of the interchange 
theorem but can be directly highlighted from the stationary character of the latter. Indeed, as the error on Z* is 
proportional to the product of the errors on the first-order change in wavefunctions, if dunk/dTK,a was known perfectly, 
a correct estimation of Z^'-^^^ should be obtained independently of the knowledge of dunw/dSfs. Putting therefore to 
zero dunk/dSp and the corresponding density changes in Eq. (41) of Ref. |l^, most of the terms cancel out and we 
recover Eq. (pT|), which evaluated for the exact dun\!_/ dr^^a must still correspond to a valid expression for Z* . 



C. First derivative of the atomic force 



By the same token as above, we can choose alternatively for 9u„k/9T„_c 
vanish in Eq. (41) of Ref. and we still obtain a valid expression for Z* 



and the associated density derivative to 



yel _ n 
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This equation corresponds to the third formulation of Z* in which it appears as the first derivative of the force on 
the atoms k, with respect to an electric field (Eq. (^)). Indeed, it is directly connected to the following expression of 
the force, deduced from the Hellmann-Fcynman theorem : 



pel _ ^0 
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S {Un\^\— —\Un]^)dk- 
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(r)][n(r)]*dr 



(14) 



Compared to Eq. ( pT| ) and Eq. (41) of Ref. ||l9|, Eq. ( p^ has the advantage that the computation of the first- 
order wavefunction derivative with respect to the electric field perturbation is the only computationally intensive step 
needed to deduce the jull set of effective charges. We note however that the implementation of Eq. (^3|), rather easy 
within a plane wave - pseudopotential approach, is not so straightforward when the basis set is dependent on the 
atomic positions, as in LAPW methods (additional Pulay terms must be introduced). 



III. BAND-BY-BAND DECOMPOSITIONS 



A. Displacement of the center of gravity of Wannier functions 

Inspired by a previous discussion by Zak Vanderbilt and King-Smith |Q emphasized that the macroscopic 
electronic polarization acquires a particular meaning when expressed in terms of localized Wannier functions. The 
periodic part of Bloch functions w„k(i") are related to the Wannier functions W„(r) through the following transfor- 
mations : 

u„k(r) = ^ E e-*-^'-^) W^{v - R) (15) 



(27r)3 

From this definition, we deduce that : 



^n(j) = ^tS^ / e^'^ '" ""k(r) dk (16) 

BZ 



d 

-Unk{ 



.„.(r) = ^ E - Wr^ir - R) (17) 
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where R runs over all real space lattice vectors. Introducing this result in Eq. (|^), we obtain 

occ „ 



(18) 



From this equation, the electronic part of the polarization is simply deduced from the position of the center of gravity 
of the electronic charge distribution, as expressed in terms of localized Wannier functions. In other words, for the 
purpose of determining the polarization, "the true quantum mechanical electronic system can he considered as an 
effective classical system of quantized point charges, located at the centers of gravity associated with the occupied 
Wannier functions in each unit cell " p4| . 

We observe that Eqs. (|l^) and (|16|) estabUsh a one-to-one correspondence between u„k and W„. As previously 
emphasized in Section III, when working within the diagonal gauge, „k becomes identified with the Kohn-Sham 
orbitals so that the associated W^^n will correspond to a single band Wannier function. Within this specific gauge, 
we can therefore isolate Pm,i3^ the contribution of band m to the f3 component of the polarization, by separating the 
different term in the sum appearing in Eq. (p^) : 



m,f3 



no 



r/3.|Wd,™(r)|2dr 



(19) 



If we take the derivative of the polarization with respect to a collective atomic displacement, Z'^ can be written 



in terms of Wannier functions as 



occ ^ 

zto.0 = Y.' J 



]dv 



(20) 



As for the polarization, this equation has a simple physical meaning. In response to an atomic displacement, the 
electronic distribution is modified and the electronic contribution to Z* can be identified from the displacement of the 
center of gravity of the occupied Wannier functions. Working within the diagonal gauge at any order of perturbation, 
we will be able to follow the change of single band Wannier functions all along the path of atomic displacements. In 
the previous expression, the contribution of band m to '^^'^ isolated : 



{Zt: 



dWm{r) 



y Wd,„(r) + (VKd,™(r)) 



]dr 



(21) 



This equation identifies the contribution from band m to the Born effective charge as fio times the change of polar- 
ization corresponding to the displacement of a point charge s on a distance equal to the displacement of the Wannier 
center of this band. Eq. (^ij) can also be estimated from finite difference by combining Eq. (^) and (p^, providing 
an easy way to decompose Z^'-^^ as soon as the Wannier functions of the system are known ||l7|. 
Alternatively, Eq. (Bl^) can also easily be evaluated in reciprocal space : 



(27r)3 



BZ 



dk 



(22) 



As Bloch and Wannier functions are related through a band-by-band transformation, the contribution from band m 
to Z* in Eq. (E2h keeps the same clear physical meaning as in Eq. (G^ 



(23) 



where Ad^ is the displacement in direction (3 of the Wannier center of band m induced by the unitary displacement of 
the sublattice of atoms k in direction a. This decomposition is strictly equivalent to what is obtained when computing 
AT'^ /3 from finite differences either in real space, using Wannier functions and Eq. (p9|), as reported by Marzari et 
al.^ or within the Berry phase approach when separating band by band contribution to Eq. (^, the reciprocal space 
equivalent of Eq. (^ . 

In practical calculations, where each band can be thought as a combinaison of well-known orbitals, the displacement 
of the Wannier center is associated to the admixture of a new orbital character to the band and must be attributed 
to dynamical changes of orbital hybridizations. As illustrated in some recent studies l9|,[l^-[l6|] , the decomposition of 
Z* appears therefore as a powerful tool for the microscopic characterisation of the bonding in solids. 
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Let us emphasize again that the previous decomposition in terms of a single band is valid only if the diagonal 
gauge was used to define the Kohn-Sham wavefunctions, hence the "d" subscript in Eq. and (|2|). The 

ground-state wavefunctions are conventionally computed within the diagonal gauge. However, in most calculations, 
the first-derivatives of these wavefunctions are computed within the parallel gauge. Within this choice, the change in 
each Bloch function will be a mixing of different Kohn-Sham orbitals when the perturbation is applied so that the 
associated change in functions Wn, defined from Eq. (jl^), will correspond to the change of a multi-band Wannier 
function. Evaluating Eq. (^) or ( |2^ ) within such a gauge, we will identify the displacement of a complex of bands 
rather than that of a single band. In practice, the first-order derivative of wavefunctions in the diagonal gauge ^^f'^U 
can be deduced from those in the parallel gauge and the ground-state wavefunctions in the diagonal gauge 

Ud,nk, by adding contributions from the subspace of the occupied bands : 



dX 



dUri 



d\ 



{Ud,nk\-9j:\^d,,nk) , . 

>, 7 ^ Ud.nk (^4j 



We note that this transformation (Eq. (|24|)) can present some problems when the denominator vanishes : this happens 
when the valence energies are degenerated. The problem can be partly bypassed by keeping a parallel transport gauge 
within the space of degenerated wavefunctions. Practically, this means that we will only be able to separate the 
contributions of disconnected set of bands. 



B. Other band- by-band decompositions 



After focusing on Eq. ([Til), we now investigate the possible band-by-band decompositions of Eq. (41) of Ref. |l9| ] 
and Eq. (|l3|). These expressions, unlike Eq. (pi]), are not written as simple sums of matrix elements, each related 
with a single band. However, individual contributions to Eq. ( |l3|) can be identified using the following decomposition 
of the density : 



n(r) 



occ 



(27r)3 Jbz^ 



s u*k(r)M„k(r)dk. 



It gives 



= 2- 



no 



occ 

^ s (u„k| 



''ext.k 



dvy,co , du„ 



for which the following decomposition is obtained, using the diagonal gauge wavefunctions 

n^ f I , c'^oxt.k , dv^co , dumi^ 



Q/3J 



'(27r)3 



S (>Ud,mk| 



BZ 



)dk 



(25) 



(26) 



(27) 



This expression corresponds to the contribution of the electrons of band m to the force induced on atom k by a 
macroscopic field £p. However, it is not equivalent to Eqs. ( ^l| ) or (p^. Indeed, for a particular band m, the 
difference between matrix elements present in Eq. (|2^) and (p7|) is (within a given gauge) : 



dv. 



\Umk.\ 



cxt.k 



dv^co . dur, 



2 



,dUr. 



,dUr. 



dkp 



(r')dr dv' 



\ I i^xc(r,r')[ 
\ f if,,(r,r')[§^(r)]*^(r')drdr' 

2 Juo ClT^a Ot/3 



(28) 



where nmk(r) is a short notation for u'^y^{r)umk{^) ■ The summation of these differences on all the bands and 
integration on the Brillouin zone gives zero, as expected. However, the band-by-band difference, Eq. ( [28[), does not 
vanish. This demonstrates that the quantity defined from Eq. ( p7| ) is independent from that of Eq. (plD and has 
therefore no specific meaning in terms of Wannier functions. Unlike Eq. (p^, Eq. (41) of Ref. |l^] is not naturally 
convertible in a sum of independent band contributions. 
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IV. NUMERICAL COMPARISON 



The previous theoretical results can now be illustrated on a numerical example. In what follows, we will consider 
the case of barium titanate (BaTiOs), a well-known ferroelectric material which is stable at high temperature in a 
cubic perovskite structure and exhibits non-trivial values of Z* jl^ . 

Our calculations have been performed within the density functional theory and the local density approximation . 
For the exchange-correlation energy, we used a polynomial parametrization [^6| of Ceperley- Alder homogeneous 
electron gas data. We adopted a planewave-pseudopotential approach. We choose highly transferable extended norm- 
conserving pseudopotentials as described in Ref. ||2^. The Ba 5s, Ba 5p, Ba 6s, Ti 3s, Ti 3p, Ti 3d, Ti 4p, O 2s, and 
02p levels have been treated has valence states. The electronic wavefunction has been expanded in plane-waves up 
to a kinetic energy cutoff of 35 hartrees. Integrals over the Brillouin zone have been replaced by sums on a 6 x 6 x 6 
mesh of special k-points. The Born effective charges have been computed in the cubic phase at the optimized lattice 
parameter of 3.94 A. They have been obtained by linear response following the scheme described in Ref. ||l9| . 

In Table H, we summarize the results obtained from independent formulations for the titanium charge (Z^,^). The 
decomposition of total Z* is provided, according to Eq. (|2|) and Eq. ( p7|) in the diagonal gauge and to Eq. (|2^) in 
the parallel gauge. We also compare our results to those reported independently by Marzari et al. [|l^ from a direct 
computation of the displacement of the center of gravity of the electronic Wannier functions. 

As expected in this class of compounds [|l6|, the total charge on the Ti atom is anomalously large (-1-7.25) and 



comparable in amplitude to the value of -1-7.16 reported independently using the Berry phase approach |29|. The 
main anomalous contribution (deviation from the nominal value of the second column) is located in the O 2p bands. 
Similarly, the oxygen charge along the Ti-0 bond is anomalously large and equal to —5.71. Both Ti and O anomalous 
charge contributions are related to each other and can be assigned to dynamical changes of hybridization between 
O 2p and Ti 3d orbitals ||l|l. This was explicitely demonstrated for a parent compound (KNbOs) by Posternak et 
al. !§. 

We observe that the global charge is equivalent independently of the approach while it is not the case for partial 
contributions coming from different isolated sets of bands. First, the band- by-band decompositions obtained within 
the diagonal and parallel gauges are not similar. This means that the unitary transform performed when changing 
the gauge strongly mix the different bands. Second, the results deduced from Eq. (^2|) and (|2^) within the diagonal 
gauge are significantly different, demonstrating that the amplitude of the quantity defined in Eq. (ES) is not negligible. 



Third, the results obtained from Eq. (22) within the diagonal gauge are comparable to those of Marzari et al. who 
explicitely computed the electronic Wannier functions and estimated Eq. ( pl| ) using a finite difference technique 
combining Eq. (^ and Eq. (^9|) . This illustrates the physical interpretation of Eq. (^2|) in terms of localized Wannier 
functions : the contributions describe the displacement of the Wannier center of each given set of bands, induced in 
response to the displacement of the Ti atom. 



V. CONCLUSIONS 



In conclusion, the Born effective charges can be computed from first-principles using different techniques and 
algorithms. The global charge is a gauge invariant quantity and is obtained independently of the approach while a 
special care is needed to separate individual contributions from separate groups of occupied bands. When using linear 
response techniques, the identification of band-by-band contributions, equivalent to those obtained within the Berry 
phase approach, requires the use of Eq. (p^, when working within the diagonal gauge. The contribution [^^'0/3]™ 
then directly related to the displacement in direction (3 of the Wannier center of band m, when displacing the sublattice 
of atoms K in direction a. The diagonal gauge condition is mandatory to identify single band contributions. The 
results obtained are conceptually and numerically different from those computed when using Eq. ( |27| ) , independently 
of the gauge choice. 
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TABLE I. Band-by-band decompositions of the Born effective charge of the Ti atom in the cubic phase of BaTiOa . The first 
hne refer to the pseudo-ion charge while the other contributions comes from the different valence electron levels. The nominal 
values expected in a purely ionic material are reported in the second column; band-by-band contributions presented in the 
three next columns were computed from linear response first-principles calculations. The last column refer to first-principles 
values deduced from the computation of Wannier functions pj] . 





Reference 




Linear response 




Wannier functions 




nominal 


Diaj 


^onal gauge 


Parallel gauge 


Diagonal Gauge 




charges 


from Eq. (p^) 


from Eq. (p7|) 


from Eq. (p2|) 


from Eq. (pl| 




+12.00 


+12.00 


+ 12.00 


+ 12.00 


+ 12.00 


Ti 3s 


-2.00 


-2.03 


+1.56 


-0.36 


-2.04 


Ti 3p 


-6.00 


-6.22 


-9.54 


-5.50 


-6.19 


Ba 5s 


0.00 


+0.05 


-0.36 


0.00 


+0.04 


O 2s 


0.00 


+0.23 


-1.56 


-0.41 


+0.20 


Ba 5p 


0.00 


+0.36 


+1.47 


+0.10 


+0.31 


O 2p 


0.00 


+2.86 


+3.68 


+1.42 


+3.01 




+4.00 


+7.25 


+7.25 


+7.25 


+7.33 
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